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Shear localization and constitutive laws in boundary region 
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We report on a numerical study of the shear flow of a simple two-dimensional model of a granular material 
under controlled normal stress between two parallel smooth, frictional walls, moving with opposite velocities 
±V. Discrete simulations, which are carried out with the contact dynamics method in dense assemblies of disks, 
reveal that, unlike rough walls made of strands of particles, smooth ones can lead to shear strain localization 
in the boundary layer. Specifically, we observe, for decreasing V , first a fluid-like regime (A), in which the 
whole granular layer is sheared, with a homogeneous strain rate except near the walls; then (B) a symmetric 
velocity profile with a solid block in the middle and strain localized near the walls and finally (C) a state with 
broken symmetry in which the shear rate is confined to one boundary layer, while the bulk of the material 
moves together with the opposite wall. Both transitions are independent of system size and occur for specific 
values of V . Transient times are discussed. We show that the first transition, between regimes A and B, can be 
deduced from constitutive laws identified for the bulk material and the boundary layer, while the second one 
could be associated with an instability in the behavior of the boundary layer. The boundary zone constitutive 
law, however, is observed to depend on the state of the bulk material nearby. 

PACS numbers: 45.70.Mg, 47.27.N-, 83.80.Fg, 83.50.Ax, 83.10.-y, 83.10.Rs 
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I. INTRODUCTION 

An active field of research over the last three decades (Ulltl, 
the rheology of dense granular flows recently benefitted from 
the introduction of robust and efficient constitutive laws. First 
identified in plane homogeneous shear flow |i3j, those laws 
were successfully applied to various flow geometries 101, such 
as inclined planes or annular shear devices lUtl, both in nu- 
merical and experimental works i^. A crucial step in the for- 
mulation of these laws is the characterization of the internal 
state of the homogeneously sheared material in steady flow 
under given normal stress by the inertial number I IHQl (see 
also Eq. ([T]i), expressing the ratio of shear time to rearrange- 
ment time, thereby regarding the material state as a general- 
ization of the quasistatic critical state, which corresponds to 
the limit of / ^ 0. Once identified in one geometry, those con- 
stitutive laws prove able to predict velocity fields and various 
flow behaviors in other situations, with no adjustable parame- 
ter jH. 

However, assuming a general bulk constitutive law to be 
available, in general, one needs to supplement it with suit- 
able boundary conditions in order to solve for velocity and 
stress fields in given flow conditions. Recent studies, mostly 
addressing bulk behavior, tended to use rough boundary sur- 
faces, both in experiments (as in Isl lioll ) and in simula- 
tions iH m [T1U13I1 . in order to induce deformation within 
the bulk material and study its rheology. Yet, in practical 
cases, such as hopper discharge flow 11411 . granular materi- 
als can be in contact with smooth walls (i.e., with asperities 



'Electronic address: [zahra. shojaaee@uni-duisburg-essen.dF 



much smaller than the particle diameter), in which case some 
slip (tangential velocity jump) is observed at the wall ifTsl - flsll . 
and the velocity components parallel to the wall can vary very 
quickly over a few grain diameters. The specific behavior of 
the layer adjacent to the wall should then be suitably charac- 
terized in terms of a boundary zone constitutive law in order 
to be able to predict the velocity and stress fields. 

In this work we use grain-level discrete numerical simula- 
tion to investigate the behavior of a model granular material 
in plane shear between smooth parallel walls^a simple setup 
which has already been observed to produce lll9l42l]l . depend- 
ing on the control parameters, several possible flow patterns, 
with either bulk shear flow, or localization of gradients at one 
or both walls. We extract a boundary layer constitutive law 
similar to the one applying to the bulk material. The stability 
of homogeneous shear profiles and the onset of localized flows 
at one or both opposite walls have been also investigated. Al- 
though Couette flow between parallel flat smooth walls is not 
an experimentally available configuration, we find it conve- 
nient as a numerical test apt to probe both bulk and boundary 
layer rheology, and their combined effects on velocity fields 
and shear localization patterns. 

The structure of the paper is as follows; Sec. HJ describes 
the model system that is simulated, and gives the definitions 
and methods used to identify and measure various physical 
quantities. In Sec. [Ill] different flow regimes are described, 
according to whether and how the velocity gradient is local- 
ized near the walls. In Sec.|IV]we derive the constitutive laws 
both in the bulk and in the boundary layer Sec.[V]applies the 
constitutive laws identified in Sec.|IV]to explain some of the 
observations of Sec. such as the occurrence of localiza- 
tion transitions or the characteristic times associated with the 
establishment of steady velocity profiles. Sec. |Vl] is a brief 
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conclusion. 



II. SYSTEM SETUP 
A. Sample, boundary conditions, control parameters 

In the contact dynamics method ll22l425ll (CD), grains are 
regarded as perfectly rigid, and the mechanical parameters rul- 
ing contact behavior are friction and restitution coefficients. 
The CD method can deal with dense as well as dilute granular 
assemblies, and successfully copes with collisions as well as 
enduring contacts, and with the formation and dissociation of 
clusters of contacting objects. We consider here a dense as- 
sembly of disks (in 2D), with interparticle friction coefficient 
/ip=0.5. As for dense frictional assemblies, the restitution co- 
efficient does not influence the constitutive laws perfectly 
inelastic collisions are considered (both normal and tangential 
restitution coefficients are set to zero (en^et—Q)). With the 
same contact properties at smooth walls (/ivi/=0.5,e„=er=0), 
slip velocities of the same order of magnitude as the shear 
velocity occur To avoid ordering phenomena, disks are poly- 
disperse, with diameters uniformly distributed between 0.8c/ 
and d. The largest diameter d is taken as the length unit 
throughout the following (d~l[L]). Similarly the mass den- 
sity of the particles is set to unity (p = l [M]/[L]^), so that 
the mass of a disk with unit diameter is m—n/A. The time 
unit is chosen such that the pressure (normal forces applied 
to the walls divided by the length of the walls) have a value 
(7„=/s./Lv=0.25[M]/[r]2, which leads to: Fy^5[M][L] / [T]^ . 
In other words, we use the following base units for length, 
mass and time: 

[L]=d, 
[M] =c/2p, 

[T] = ^Sd^p/Fy. 

We consider simple shear flow within rectangular cells with 
periodic boundary conditions in the flow direction (parallel 
to the X axis in Fig. [U. Gravity is absent throughout all our 
simulations. The top and bottom walls bounding the cell are 
geometrically smooth, but their contacts with the grains are 
frictional, with a friction coefficient /Zw set to 0.5. They move 
with constant and opposite velocities {±V) along direction x. 
They are both subjected to inwards oriented constant forces Fy 
normal to their surface, so that in steady state a constant nor- 
mal stress Uyy is transmitted to the sample. The wall motion in 
the normal direction is ruled by Newton's law, involving the 
wall mass, equal to 50, thereby causing the system height L, 
to vary in time. In steady state shear flow, Ly fluctuates about 
its average value. 

Results from different samples of various sizes are pre- 
sented below. System sizes and simulation parameters are 
listed in Tab. U 

As in Refs. IH-Ht], the dimensionless inertial number, /, de- 




FIG. 1 : (color online) A polydisperse system of hard frictional disks 
in planar shear geometry with periodic boundary conditions in x di- 
rection. A prescribed normal force Fy to the confining walls, deter- 
mines the constant external pressure of the system. The walls move 
with the same constant velocity V in opposite directions. The width 
of the lines connecting the centers of contacting particles represents 
the magnitude of the normal contact forces above a threshold. 
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511 
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620 


20000 
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1023 


40 


20 
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0.03-30.00 


2500 


10000 
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1023 


40 


20 


0.0625 


0.03-30.00 


9900 


10000 
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3199 


50 


50 


0.25 


0.01-30.00 


4000 


8000 
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2047 


80 


20 


0.25 


0.01-20.00 


10000 


4000-12000 
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3071 


120 


20 


0.25 


0.01-35.00 


22000 


6000 


7 


5119 


200 


20 


0.25 


0.01-30.00 


64000 


13000 



TABLE I: Parameters used in the simulations, n is the number of 
disks in the sample. Tjs denotes the characteristic time to approach 
steady state according to Eq. ( 117b . T^ira is the total time simulated in 
each run. 



fined as a reduced form of shear rate j: 




is used to characterize the state of the granular material in 
steady shear flow. In contrast to previous studies iHQl, the 
shear is not homogeneous in the present case (because of wall 
slip and of stronger gradients near the walls), and in general 
2V 

y is different from — . Thus, the shear rate has to be mea- 

Ly 

sured locally. We focus in the present study on shear localiza- 
tion at smooth walls and try to deduce constitutive laws in the 
boundary layer, associating the boundary layer behavior not 
only with wall slip, but also with the material behavior in a 
layer adjacent to the wall, the internal state of which might be 
affected by that of the bulk material. 



B. System preparation 

To preserve the symmetry of the top and the bottom walls, 
the system is horizontally filled. While distance Ly between 
the walls is kept fixed, a third, vertical wall is introduced, on 
which the grains (which are temporarily rendered frictionless) 
settle in response to a "gravity" force field parallel to the x 
axis. Then the force field is switched off, and the free sur- 
face of the material is smoothened and compressed by a pis- 
ton transmitting a^x = 0-25 (the same value as o\,y imposed in 
shear flow), until equilibrium is approached. System width Lx 
is determined at this stage. Then the vertical wall and the pis- 
ton are removed, the friction coefficients are attributed their 
final values /Xp and pLw and periodic boundary conditions in 
the X direction are enforced. With constant L, and variable 
Ly, the shearing starts with velocities ±.V for the walls and an 
initial linear velocity profile within the granular layer. 



C. Measured quantities 

Before presenting the results, we first explain the method 
used to measure the effective friction coefficient, the velocity 
profiles and the inertial number. For a system in steady state, 
assuming a uniform stress tensor in the whole system, a com- 
mon method to calculate the effective friction coefficient is to 
average the total tangential and normal forces acting on the 
walls over time and then calculate their ratio. Another way 
to calculate the effective friction coefficient is to consider the 
components of the stress tensor with its contact, kinetic and 
rotational contributions inside the system [5i,26i]- The stress 
in our system is dominated by contact contributions. Let 0^ 
denote the total contact stress tensor calculated for each parti- 
cle / with aieaAi—ndf /A: 

< = J,L'^ij®r-Ij. (2) 

The summation runs over all particles j having a contact with 
particle /. F^j is the corresponding contact force and rjj de- 
notes the vector pointing from the center of particle / to its 
contact point with particle j. We used both methods, but find- 
ing no significant difference, we present in all our correspond- 
ing graphs the effective friction coefficient (/ieff) measured in 
the interior of the system considering all terms of the stress 
tensor, although the contact contribution dominates. 

Our calculation of the velocity profile accounts for parti- 
cle rotations, which contribute to the local velocities averaged 
in stripes of thickness Ay=l along the flow direction, as fol- 
lows. To each horizontal stripe centered at y~y', we attribute 
a velocity by averaging the contributions of all the particles it 
contains (partly or completely) IztIi : 

J {'Oix + COiriy)dS 
M/)- ' y.,. ■ (3) 
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FIG. 2: (color online) Transient to steady state for V = 0.70 in a 
system with Lx = 50 and Ly = 50. 

Sj denotes the surface fraction of particle / within the stripe, 
Vix its center of mass velocity in x direction, CO, its angular 
velocity and r,y is the vertical distance between the center of 
mass of the particle and a differential stripe of vertical position 
y and surface dS within surface Sj. The velocity profiles pre- 
sented here are also averaged over time intervals of Af=80. 
Those time intervals follow each other directly without any 

gap- 
In the calculation of the profiles of stress tensor, each parti- 
cle contributes to each stripe in proportion to the surface area 
contained in the stripe. This corresponds to the scheme used 
in ll27ll and is slightly different from the coarse graining re- 
viewed in ll28ll in the sense that it is highly anisotropic (with a 
coarse graining scale of Ly x 1 ) and does not incorporate the 
(stress free) regions beyond the walls. One other method is to 
split the contact contributions proportionally to their branch 
vector length within each stripe. One may also cut through 
the particles and add up the contact forces of all cut branch 
vectors. All three different methods lead to the same results 
in our simulations. 



III. VELOCITY PROFILES AND STRAIN 
LOCALIZATION 

A. Steady state 

A system sheared with a certain constant velocity under 
prescribed normal stress is expected to reach a steady state 
after a transient. For instance, in a system of size L^=50 and 
Ly=50 with a large shear velocity, V=QJ, the steady state is 
reached after a shear distance of about A c± 420, correspond- 
ing to a shear strain of 7 ~ 8 (Fig. |2]). The shear distance 
is calculated by multiplying the total shear velocity (2V) by 
time. Due to slip at the smooth walls and because of non- 
homogeneous flow, the values attributed to the shear distance 
and the shear strain overestimate the real values in the bulk 
material. Transient times before steady state will be estimated 
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FIG. 3: (color online) (a) Center of mass velocity versus shear strain 
for y = 0.70 in a system with Lj=50 and Lj,=50. The dashed red 
lines represent the velocity of the top and bottom walls, (b) His- 
togram of center of mass velocity (accumulated over a long time and 
over different simulated systems). 



in Sec. IV Al based on the constitutive laws. 

In the steady state, the center of mass velocity in the system 
of Fig. [3]is not conserved, because of the constant velocity of 
the walls, and fluctuates about its vanishing average with an 
amplitude amounting to about 10% of velocity V. This rela- 
tively high level of fluctuation, which is expected to regress in 
larger samples, reflects granular agitation within the sheared 
layer at large Y . The fluctuations in height Ly and solid frac- 
tion V (measured in the whole system) amount to only about 
1% of the average after a short transient (Fig. |4]i. The ini- 
tial sharp drop of v, for very small shear strains, is due to the 
combined effects of shear flow onset, friction activation and 
change of boudary conditions on the configuration prepared 
as described in Sec. Ill Bl 

In the steady state the profiles of the effective friction co- 
efficient stay almost uniform throughout the system (Fig. |5]l, 
but fluctuate in time, which is a direct consequent of shear- 
ing with constant velocity and consequent fluctuations in the 
center of mass velocity. 



B. Shear regimes and strain localization 

Fig.|6]displays the time evolution of the velocity profiles of 
a system of initial height Ly = 1 20 sheared with different veloc- 
ities (a) y =2.0, (b) y=0.2 and (c) V =0.03. Those three cases 
are characteristics of three different regimes observed in dif- 
ferent intervals as velocity Y decreases. At large velocity V, 
as for y=2.0 (panel (a) in Fig.|6]l, the velocity profile adopts 
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FIG. 4: Solid fraction v versus shear strain at V=0.70 in a system 
with L;t=50 and Ly=50. 




FIG. 5: (color online) Profiles of the measured effective friction co- 
efficient at different times in steady state for V=0.70 (Lv=50 and 
L,,=50). 



(after a transient) a symmetric, almost linear shape with only 
small fluctuations in time. The shear rate is somewhat larger 
near the walls than in the bulk, but the latter region is ho- 
mogeneously sheared. This is the fast or homogeneous shear 
regime (regime A in the sequel). For intermediate velocities, 
such as y=0.2 (panel (b) in Fig.|6]l the shear rate is strongly 
localized at the walls, while, ten grain diameters away from 
the walls, the material is hardly sheared at all. While the 
profile shape is essentially stable, its position on the veloc- 
ity axis fluctuates notably: the bulk material behaves like a 
solid block, but its velocity exhibits large fluctuations. To this 
situation we shall refer as the intermediate or two-shear band 
regime (regime B). Finally, for a low enough shear velocity, 
as for y=0.03 (panel (c) in Fig.|6ll the profiles fluctuate very 
strongly and the top/bottom symmetry is broken. The shear 
strain strongly localizes at one wall, while the rest of the sys- 
tem, the bulk region and the opposite wall, moves like one sin- 
gle solid object lfT9l - l2lll . This is the slow shear or one-shear 
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FIG. 6: (color online) Velocity profiles at different shear distances 
for three different values of V, as indicated, in sample with height 
Lv=120 (System 6 in Tab.|Ill. 



hand regime (regime C). Localization occasionally switches 
to the other wall, with a transition time that strongly depends 
on the system size and on the shear velocity (as we shall dis- 
cuss in Sec. IV Al ). 

In regime A the sheared layer behaves similarly to the ob- 
servations reported by da Cruz et al. iH, in a numerical study 
of steady uniform shear flow of a granular material between 
rough walls. However, with rough walls the homogeneous 
shear regime persists down to very low velocities, in spite of 
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FIG. 7: (color online) Profiles of velocity and effective friction coef- 
ficient (inset) in steady state and in the transient states for V=0.08 in 
a system with Lv=50 and Lv=50. 



increasing fluctuations. The smooth walls of our system, al- 
lowing for slip and rotation at the walls, are responsible for 
the more complex behavior lfl9l - l2lll . 

In order to be able to observe the three regimes, sample 
height Ly should be large enough. In smaller systems {Ly < 80) 
the effects of the boundary layers on the central region are 
strong enough to preclude the observation of a clearly devel- 
oped intermediate regime. Sheared granular layers of smaller 
thickness most often exhibit a direct transition from regime A 
to regime C on decreasing velocity Y . 

Our system size analysis shows a discontinuous transition 
from regime B to regime C, at Vbc— 0.10 and a continuous 
transition between regimes A and B completed at y^B~0.50 
lfT9l - l2lll . Vbc and V^b are system size independent. 

Upon reducing the shear velocity in the intermediate shear 
regime towards Ybc larger and larger fluctuations in the ve- 
locity fields are observed, involving increasingly long corre- 
lation times. Slightly above Y^q the approach to a steady state 
becomes problematic, even after the largest simulated shear 
strain (or wall displacement) intervals. Then below Ybc the 
width of the distribution of the bulk region velocities reaches 
its maximum value, TV , and the velocity profile stays for 
longer and longer time intervals in the localized state with one 
shear band at a wall (regime C). Such localized profiles can be 
regarded as quasi-steady states - as switches from one wall to 
the opposite one, ever rarer at lower velocities, sometimes oc- 
cur. The lifetime of these one-shear band asymmetric steady 
shear profiles also increases with system height Ly, similar 
to ergodic time in magnetic systems I19i, i20il . These quasi- 
steady states also exhibit uniform stress profiles, contrary to 
the nonuniform ones in the transient states, as the localization 
pattern is switching to the other side (Fig.|2l). 

Fig.[8](a) is a plot of center of mass velocity in the flow di- 
rection versus time in regime C. Most of the time, it is slightly 
fluctuating about the value of either one of the velocities of the 
walls, ±y, as also illustrated in the histogram plot, Fig.[8](b), 
for which values were accumulated over a long time and over 
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FIG. 8: (color online) (a) Center of mass velocity fluctuations in 
steady state for K=0.05 in a system with L^=20 and L,.=20. The 
dashed red lines represent the velocity of the top and bottom walls. 
The transition time (magnified in the inset) is measured at both ends 
of direct transitions from one wall to the other, between the round 
dots, (b) Histogram of center of mass velocities. 
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FIG. 9: (color online) Slip velocity u^"P (averaged over and 
U^''') measured as a function of shear velocity (systems specified in 
TaklB. 



general tendency is the same, slight deviations are observable. 
The apparent change of slope of the graph near V=l could 
be associated with larger strain rates and inertial numbers in 
boundary layers, gradually approaching a collisional regime 
(see Sec. lIVCl about the coordination number). 



IV. CONSTITUTIVE LAWS 



different simulated systems. Transition times as the shear 
band switches directly from one wall to the other are mea- 
sured as indicated. Those times are recorded to be discussed 
in Sec. [VA2l 

C. Slip velocity 

The slip at smooth walls is a characteristic feature of the 
boundary region behavior. 

To evaluate the slip velocity at the walls one needs to calcu- 
late the average of the surface velocity of particles in contact 
with the walls at their contact point. The slip velocity in this 
work is defined as the absolute value of the difference between 
the wall velocity and the average particle surface velocity at 
the corresponding wall, u^''^ at the bottom, respectively v^^'^ 
at the top wall. To this end all particles in contact with the 
walls over the whole simulation time in steady state should be 
considered, and contribute 

Uo''P = y + (u,, + o),T,),v, (4) 
V^;'P = y-(u,-,-co;r;),.,, (5) 

where u,v is the x component of the center of mass velocity of 
particle / of radius r, with angular velocity O),-. 

Our observations show that the slip velocity in a certain 
shear velocity interval 0.2 <y< 1.0 does not depend on the 
system size (Fig. |9]l. For larger shear velocities, though the 



Constitutive laws were previously studied, in similar model 
materials, in homogeneous shear flow iHIZ^IlSl- Their sensi- 
tivity to material parameters (restitution coefficients, friction 
coefficients and, possibly, finite contact stiffness) is reported, 
e.g. in Isill . In our system we separate the boundary re- 
gions near both walls, from the central one (or bulk region). 
Unless otherwise specified, the boundary regions have thick- 
ness /i = 10. Near the walls, the internal state of the granular 
material is different, and we seek separate constitutive laws 
for the boundary layers and for the bulk material. While the 
bulk material is expected to abide by constitutive laws that 
apply locally, and should be the same as the ones identified 
in other geometries or with other boundary conditions IHIst], 
the boundary constitutive law is expected to relate stresses to 
the global velocity variation across the layer adjacent to the 
wall. In a continuum description suitable for large scale prob- 
lems, this will reduce to relating stresses to tangential velocity 
discontinuity. 



A. Constitutive laws in the bulk region 

1. Friction law 

The steady state values of the inertial number (/bulk) and 
that of the effective friction coefficient figff are measured, as 
averages over time and over coordinate y within the interval 
li < y < Ly — h. jXeff is plotted as a function of /bulk for all 
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FIG. 10: (color online) /igff as a function of inerdal number in the 
bulk region for different system sizes (see Tab.|lJ. The fit function is 
calculated according to Eq.|6]for ^0=0.25. The error bars are much 
smaller than the symbols. The inset is a semilogarithmic plot of the 
same data. 
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FIG. 1 1 : (color online) Influence of h on fl^fi as a function of inertial 
number in the bulk region (data from systems 5 and 7 in Tab.|ll(. The 
dashed red line represents the critical friction coefficient /io=0.25. 



FIG. 12: (color online) v as a function of inertial number in the bulk 
region. The error bars plotted are much smaller than the symbols 
(systems specified in Tab.|D. 



terial to be continuously sheared, the ambient noise level, due 
to the proximity of the sheared boundary layer, entails slow 
rearrangements that produce macroscopic shear ll32ll . Upon 
increasing h the central bulk region excludes the outer zone 
that is affected by this creep effect. The critical friction coef- 
ficient, from Fig.[TTl is /Xo=0.25 (below which the data points 
are sensitive to the value of h), which is consistent with the 
results of the literature ^ [st] . 

Fitting /ieff — Ik) with a power law function, as in l29l[30l] 



Meff-MO=A-/j 



bulki 



(6) 



the following coefficient values yield good results (see 
Fig.IIOli: 

;Uo = 0.24 ±0.01, 
A = 0.92 ±0.05, 
5 = 0.80 ±0.05. 



2. Dilatancy law 



different system sizes in Fig. [TO] showing data collapse for 
different sample sizes. 

The apparent influence of the choice of h on the measured 
effective friction coefficient and inertial number in the bulk 
region is presented in Fig. [TT]for two different system sizes 
and for two different h values. 

We observe that some data points with finite values of /bulk 
(4uik > 10^^) are shifted to much smaller values of /bulk upon 
increasing h: compare the open and full symbols in Fig. [TT| 
This effect is apparent in regimes B and C. It is due to the 
creep phenomenon (as was also observed in the annular shear 
cell in IH), which causes some amount of shearing at the 
edges of the bulk region, adjacent to the boundary layer, al- 
though the effective friction coefficient is below the critical 
value. Although the local shear stress is too small for the ma- 



We now focus on the variation of solid fraction v as a func- 
tion of inertial number within the bulk region, v is averaged 
over time, once a steady state is achieved, within the cen- 
tral region, h < y < Ly — h. Function Vbuik(4uik) is plotted 
in Fig. [12] for different system sizes, leading once again to a 
good data collapse. A linear fit for all data sets in the interval 
0.03 < /bulk < 0.20 gives: 

Vbuik = 0.81 -0.30 -/bulk, (7) 
which is consistent with the linear fit in IHIsl]. 

B. Constitutive laws in the boundary layer 

In order to characterize the state of the boundary layer of 
width h adjacent to the wall (recall h = 10 by default), we use 
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FIG. 13: (color online) /Xeff a function of inertial number in 
the boundary layer. The error bars plotted are much smaller 
than the symbols. As /boundaiy >4iilk (shear localization at smooth 
walls) Ateff (boundary) lies always beneath Ateff(4ulk)- The presented 
/leif(/bulk) curve is based on the same data set as Fig. 1101 



a local inertial number /boundaiy. defined as follows: 



top/bottom^ Al^„ M^top/bottom 
boundary W _ \ 



with 



(8) 



(9) 



Au'°P = y-u,(Lv-/z), 

Au''°"°'^ = -u,-(/!) + y. 

1. Friction law 



Fig- mis a plot of /ieff as a function of the inertial number 
boundary in the boundary layer for all different system sizes. 

In steady state the value of /igff in the boundary layer has to 
be equal to the averaged one in the bulk. The observed shear 
increase (in regime A) or localization (in regimes B and C) 
near the smooth walls entails larger values of inertial numbers 
in the boundary region. An equal value of /igff in the bulk 
and in the boundary zone then requires that the graph of func- 
tion /ieff (boundary) is below its bulk Counterpart in the inertial 
number interval measured. 

In Sec. lIV Al we have seen that the friction law can be iden- 
tified in the bulk independently of h (see Fig.fTTTi. as an intrin- 
sic constitutive law. According to the definition of /boundary 
in Eqs. dHJ and (|9]l any constitutive relation involving /boundary 
should trivially depend on h. In shear regimes B and C, there 
is no shearing in the bulk region, and consequently Au in the 
numerator of Eq. ([8]l does not change with h. On multiplying 
the measured /boundary with the corresponding value of h, we 
thus expect the data points belonging to shear regimes B and 
C to coincide (Fig.[T4]i. In regime A, in contrast, the existence 
of shear in the bulk region leads to an apparent h dependence 
of the measured Au. Accordingly, after multiplying /boundary 



with /i, the curves do not merge. The critical effective friction 
coefficient at which the deviation of the curves begins corre- 
sponds to /io=0.25 (the dashed horizontal line in Fig.O, in 
agreement with the results in Sec. II V A II This makes it more 
difficult to identify a constitutive law for the boundary layer, 
when the bulk region is sheared in regime A. 
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FIG. 14: (color online) /ieff versus h x /boundary on linear (top panel) 
and semi-logarithmic (bottom panel) plots. The dashed horizontal 
line indicates the critical state value /i^ff = /io=0.25. 

The behavior of /ieff shown in Fig.[T4]is apparently anoma- 
lous in two respects: (;) the Au dependence of /igff does not 
seem to follow a single curve (suggesting /ieff depends on 
other state parameters than the velocity variation across the 
boundary zone); (//) /ieff is a decreasing function of /boundary 
for the first data points, as /; x /boundary < 0-2. In Fig. 1151(a). 
we take a closer look at the low /boundary data points, which 
bear number labels 1 to 6 in the order of increasing shear ve- 
locity Y . The transition from regime C (one shear band) to 
regime B (two shear bands) occurs between points 4 and 5, 
whence a decrease in /boundary- as the velocity change across 
the sheared boundary layers changes from TV to merely Y . In 
an attempt to identify one possible other variable influencing 
boundary layer friction, the symbols on Fig.[T5](a) also encode 
the value of the bulk density. We note then that points 4 and 6, 
which have different friction levels, although approximately 
the same /boundary- correspond to different bulk densities. The 
constitutive laws in the boundary layer might thus depend on 
parameter Vbuik in addition to /boundary- 

As to issue (//), the decrease of /ieff before the zig-zag pat- 
tern on the curve of Fig.[T5](a) (data points 1 to 3) is associ- 
ated to an increase in the boundary layer density with /boundary - 
This is not the case in all of the systems and these features 
strongly depend on the preparation and the initial packing den- 
sity (compaction in the absence of friction). Independent of 
whether /igff in regime C increases or not as /boundary increases, 
/Ieff always displays a decreasing tendency as Vboundary in- 
creases, just like /ieff and v vary in opposite directions in bulk 
systems under controlled normal stress, as shown in Ref. jstl, 
or as expressed by Eqs. (|6]l and (|7]i (panel (b) in Fig.fTSll. The 
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FIG. 15: (color online) (a) ^eff as a function of /boundary (data from 
system 5).The full symbols belong to the states with bulk densities 
larger than the critical value Vc=0.81 (see Eq. (O). The full dia- 
monds have a density between 0.81 and 0.82, full squares have a 
density between 0.82 and 0.83 and full circles have a density larger 
than 0.83. The inset represents the zig-zag with some more data 
points, which are absent in the master graph for the sake of clarity, 
(b) jlgff as a function of Vboundary- The error bars are smaller than the 
symbols. 



lack of a perfect collapse of the data points around the de- 
creasing linear fit of Fig. [15] (b) shows however that the state 
of the boundary layer in slowly sheared systems does not de- 
pend on a single local variable, but is influenced by the state 
of the neighboring bulk material, as remarked above. 



2. Dilatancy law 

After averaging the profiles of solid fraction and inertial 
number over the whole simulation time in steady state in the 
boundary region, Vboundary (4oundary) graphs ai-e then plotted in 
Fig.[T6](a) for different system sizes. In Fig.[T6](b). Vbuik(-4uik) 

and Vboundary (boundary) are Compared for all data sets. Vboundary 

and Vbuik drop proportionally with increasing /bulk (shear ve- 
locity) until /bulk— 0.08 (in shear regime A). Afterwards, the 
drop in Vboundary is much Steeper (Fig. [16] (c)). 



FIG. 16: (color online) (a) v as a function of inertial number in the 
boundary layers (systems specified in Tab.|Ill. (b) Vboundary (boundary) 
compared to Vbuik(4ulk)- (c) The ratio between Vboundary and Vbuik as 
a function of /bulk- 



C. Coordination number 

The coordination number (average number of contacts per 
grain) is a quantitative measure of the status of the contact 
network. Fig.[T7]shows the measured coordination number in 
the bulk and in the boundary layers as a function of inertial 
numbers in these two regions. The data are collected from 
different systems in Tab. [I] 

The bulk coordination number, Zbuik is fitted with the power 
law function Zbuik = 2.70 — 2.76 • /bulk- boundary re- 

gion coordination number, Zboundai^, follows a slightly dif- 
ferent dependency on /boundary; which becomes noticeable for 
boundary ^0.1, which Corresponds to y~l (see Sec lIIICI l. It 
drops to smaller values, as /boundary reaches larger values, 
above 1 . The decrease of coordination number as a function of 
inertial number is compatible with the observations of lid, in 
which some effect of restitution coefficient on Z was however 
reported. The finite softness of the particles is also known to 
affect coordination numbers ifllll much more than the rheo- 
logical laws. It is only for configurations extremely close to 
equilibrium that coordination numbers are observed to exceed 
the minimum value 3 for stable packings of frictionless disks 
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FIG. 17: (color online) Coordination number as a function of inertial 
number in the bulk and boundary regions. The error bars are much 
smaller than the symbols. 



(excluding rattlers), and even in quite slow flows this condi- 
tion is not fulfilled. 



V. APPLICATIONS 

We now exploit the constitutive relations and other observa- 
tions reported in the previous sections to try and deduce some 
features of the global behavior of granular samples sheared 
between smooth walls. 



A. Transient time 



1. Transient to steady state in regime A 



The bulk friction law of Sec. llV ATl can be used to evaluate 
the time for a system to reach a uniform shear rate in regime A, 
if we assume constant and uniform solid fraction v and normal 
stress Oyy, and velocities parallel to the walls at all times. We 
write down the following momentum balance equation: 



dt 



dy 



(10) 



looking for the steady solution: Vx = jy- Assuming constant 
p, V and Oyy we can write: 



dv^ d 
which leads by derivation to: 



[Meff(7)] (yyy 



(11) 



(12) 



Separating the shear rate field into a uniform part 70 and a y- 
dependent increment A7, and assuming as an approximation 



just a linear dependency of /ieff on 7, we can rewrite Eq. ( fT2b 
as follows: 



dAj ^Meff . . 
which is a diffusion equation with diffusion coefficient 

dHe« CTyy 



D 



df pv' 



(13) 



(14) 



The characteristic time to establish the steady state profile 
(uniform 7 over the whole sample height Ly) is then: 



7ss = -^. 



(15) 



A linear fit of function /ieff(4urk) (see Fig. [TOb in interval 
(0.03 < /bulk < 0.20) is: 

Meff = 0.27 +1.16 -/bulk. (16) 
According to Eqs. O, (O, (Ol and this leads to: 

rss^l.56L2. (17) 

The estimated values Tss for different system sizes is listed 
in Tab.|Il As T^s grows like L^, very long simulation runs be- 
come necessary to achieve steady states in tall (large Ly) sam- 
ples, and some unstable, but rather persistent, distributions of 
shear rate can be observed Isol Q. Our data for Lv=120 
and Lv=200 may still pertain to slowly evolving profiles, even 
though the constitutive law can be measured in approximately 
homogeneous regions of the sheared layer over time intervals 
in which profile changes are negligible. 

2. Transition from one wall to the other in regime C 



As stated in Sec. 1111 B I in regime C the asymmetric veloc- 
ity profiles can be regarded as steady states and the switch- 
ing stages in which the shear band changes sides are tran- 
sient states in which the shear stress is not uniform through 
the granular layer We now try to estimate the characteristic 
time for such transitions. This estimation does not rely on 
a specific model for the triggering mechanism of the transi- 
tion. It is based on the simple idea that the transition takes 
place when the solid block is accelerated due to a shear stress 
difference between the top and the bottom boundary zones. 
Taking the whole bulk region as a block of mass M moving 
with the velocity of the top wall V, a transition to velocity —V 
with acceleration A will take: 



_2V 

^transition — ~;~ i 
A 

in which the acceleration A is equal to: 
(err 



A = 



M 



(18) 



(19) 
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FIG. 18: (color online) Transition time divided by the shear velocity 
as a function of system height. Empty circles correspond to measured 
times, while the (red) dashed line plots estimated ones, using ( |20| l. 



FIG. 19: (color online) The critical /boundary which corresponds to 
/i()=0.250 (dashed red line) and determines the the critical velocity 
Vab for the transition from regime A to regime B. 



Substituting M=p\L„Ly and a'^ - CJ^f """=Ajti(j,.,: with 

.top _ bottom , 



one gets: 



2pvVLy 

AflGyy 



(20) 



Accordingly, the transition time increases proportionally to 
the shear velocity and to system height L,,. Using v ~ 0.84, 
Gyy = 0.25 and taking Aji ~ 0.05 as a plausible value in shear 
regime C (see Figs. l7land[T4li we calculate as a func- 

tion of system height L,.. In Fig.[T8]these estimated times are 
compared to transition times that are measured as explained 
in the caption of Fig. [8] 

Admittedly, one does not observe only direct, sharp tran- 
sitions in which localization changes from one wall to the 
opposite one. Some transient states are more uncertain and 
fluctuating, and the system occasionally returns to a localized 
state on the same wall after some velocity gradient has tem- 
porarily propagated within the central region. The data points 
of Fig. [TS] correspond to the well-defined transitions, which 
become less frequent with increasing system height. Thus a 
unique data point was recorded for systems with Ly = 120 and 
Ly = 160. The comparison between estimated and measured 
transition times is encouraging, although the value of A/i in 
(I20I 1 is of course merely indicative (it is likely to vary during 
the transition), and the origin of such asymmetries between 
walls is not clear. 



B. Transition velocity Vab 

/io=0.25 from the power law fit in Eq. ^ corresponds to 
the minimal value of the bulk effective friction coefficient, 
the critical value below which the granular material cannot 
be continuously sheared (except for local creep effects in the 
immediate vicinity of an agitated layer). 

Fig.[T9]gives the value of the inertial number in the bound- 
ary region, such that the boundary friction coefficient matches 



Alo=0.25: 



Mo = 0.25 /boundaiy=0.086 ±0.005. 



(21) 



Thus for /boundary ^S0 086 we expect no shearing in the bulk. 
According to Eqs. ^ and ^ this results in y=0.485 ±0.028, 
in very good agreement with our observations reported in 
Sec.|inB](V:^B~0.50). 

The explanation of the transition from regime A to regime B 
is simple: the boundary layer, with a smooth, frictional wall, 
has a lower shear strength (as expressed by a friction coef- 
ficient) than the bulk material. Thus for uniform values of 
stresses (7„. and o;„. in the sample, such that their ratio C7,-,./ Gyy 
is comprised between the static friction coefficient of the bulk 
material and that of the boundary layers, shear flow is confined 
to the latter. 



C. Transition to regime C at velocity Vbc 

Although it is not systematically observed, and is likely to 
depend on the bulk density, the decreasing trend of /ieff in the 
boundary layer as a function of Av or of /boundaiy, as appar- 
ent in Figs. [T4land[T5](a). provides a tempting explanation to 
the transition from regime B to regime C. Assuming /igff for 
given, constant Gyy, to vary in the boundary layers as 



Meif = Mo - witha>0, 



(22) 



one may straightforwardly show that the symmetric solution 
with Au = ±y, and solid bulk velocity = 0, is unstable. A 
simple calculation similar to the one of Sec. lV A2l shows that 
velocity Vs, if it differs from zero by a small quantity 8Vs at 
r = 0, will grow exponentially. 



X),(t) = Sujexp — 



(23) 



12 



until it reaches ±V, with the sign of the initial perturba- 
tion 8Vs- Transition velocity Vbc would then be associ- 
ated to a range of velocity differences Av across the bound- 
ary layer with softening behavior (i.e., decreasing function 

Meff (^boundary))- 

In view of Fig. [15] (a), where the BC transition takes place 
at point 4, this seems plausible, as the slope of function 
Meff (boundary) appears to Vanish towards this point. 

VI. CONCLUSION 

In this work we have investigated shear localization at 
smooth frictional walls in a dense sheared layer of a model 
granular material. The slip at the walls induces inho- 
mogeneities in the system leading to three different shear 
regimes. As the wall velocity is reduced from large values, 
two transitions successively occur, in which shear deforma- 
tion localizes, first symmetrically near opposite walls, and 
then at a single wall. Measuring stress tensor, inertial number 
and solid fraction locally in the whole system, the constitu- 
tive laws have been identified in the bulk (for which our re- 
sults agree with the published literature) and in the boundary 
layer Those constitutive laws, supplemented by an elemen- 
tary stability analysis, allow us to predict the occurrence of 



both transitions, as well as characteristic transient times. The 
consistence of the derived constitutive laws for the bulk rhe- 
ology with those in previous contributions Jst] using the MD 
method, confirm that the rheology is the same for CD and for 
MD in the limit of large contact stiffness. 

Additional numerical work should be carried out in order 
to assess the dependence of the boundary layer constitutive 
law on the state of the adjacent bulk material with full gener- 
ality. The application of similar constitutive laws for smooth 
boundaries should be attempted in a variety of flow configura- 
tions; inclined planes, vertical chutes, circular cells. Finally, 
the success of the simple type of stability analysis carried out 
in the present work calls for more accurate, full-fledged ap- 
proaches in which couplings of shear stress and deformation 
with the density field would be taken into account. 
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